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Abstract. We present new results on dynamical instabilities in rapidly rotating neutron- 
stars. In particular; using numerical simulations in full General Relativity, we analyse the 
effects that the stellar compactness has on the threshold for the onset of the dynamical bar- 
mode instability, as well as on the appearance of other dynamical instabilities. By using an 
extrapolation technique developed and tested in our previous study we expHcitly determine 
the threshold for a wide range of compactnesses using four sequences of models of constant 
baryonic mass comprising a total of 59 stellar models. Our calculation of the threshold is 
in good agreement with the Newtonian prediction and improves the previous post-Newtonian 
estimates. In addition, we find that for stars with sufficiently large mass and compactness, 
the m=3 deformation is the fastest growing one. For all of the models considered, the non- 
axisymmetric instability is suppressed on a dynamical timescale with an m = 1 deformation 
dominating the final stages of the instability. These results, together with those presented 
in [T], suggest that an m = l deformation represents a general and late-time feature of non- 
axisymmetric dynamical instabilities both in full General Relativity and in Newtonian gravity. 



PACS numbers: 04.25.Dm, 04.30.Db, 04.40.Dg, 95.30.Lz, 95.30.Sf 97.60.Jd 

1. Introduction 

Non-axisymmetric deformations of rapidly rotating self-gravitating bodies are rather generic 
phenomena in nature and could appear in a variety of astrophysical scenarios like stellar core 
collapse 111, accretion-induced collapse of white dwarfs |3|, or the merger of two neutron 
stars H. Over the years a considerable amount of work has been devoted to the search of 
unstable deformations that, starting from a quasi-axisymmetric stellar configuration, would 
lead to the formation of highly deformed rotating massive objects (see ref. HI for a detailed 
list of references). One of the main reasons behind this interest is that such deformations 
would lead to the intense emission of high-frequency gravitational waves {i.e., in the kHz 
range) which is potentially detectable by ground-based detectors such as LIGO, GEO, Virgo 
or the planned resonant detector such as DUAL IS). 

Despite such extensive studies, various questions about the dynamics of the non- 
axisymmetric deformation of rapidly rotating self-gravitating bodies are not yet completely 
clarified. Among the most important questions that have been addressed only rather recently 
it is worth to recall the following ones: i) How long do these deformation survive once 
they reach their maximum amplitude? ii) How large is the energy emitted in gravitational 
waves? Hi) Which physical phenomena determine the shortest damping timescale and impress 
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a signature on the emitted signal? iv) What is the effect of the stellar compactness M / R^, 
where M and Re are the stellar mass and the proper equatorial radius, respectively, on the 
dynamics of the instability and on the threshold for its onset? 

While points i)-iii} were first addressed in ref. LIJ (hereafter paper I), here we concentrate 
on providing an answer to question iv) supplying new information on the general properties 
of the dynamical instability in a very large class of stellar models that are characterized by 
differential rotation and by high compactness and that are members of four sequences of 
models with constant baryonic mass. 

The main result obtained is that any non-axisymmetric deformation that develops in 
our models is damped over a dynamical timescale, through pure inviscid hydrodynamical 
nonlinear phenomena. Moreover, for all the models that develop dynamical instabilities, the 
m=l deformation eventually becomes the dominant one irrespective of whether the models 
are above or below the threshold for the development of the bar-mode instability (see also 
paper I). This evidence is consistent with the simulations performed in Newtonian gravity 
by Ou and Tohline for stars with a very-low (3 (where (3 = T/|M^| is the ratio between 
the rotational kinetic energy T and the gravitational binding energy W), thus suggesting that 
this may be a general feature of this type of dynamical instabilities. In addition we show 
that, in a region of high stellar compactness, other instabilities, such as one having an m=3 
deformation, can develop. Finally, adopting an extrapolation technique developed and tested 
in paper I, we determine the threshold pc for the onset of the bar-mode instability for all of the 
sequences considered, thus determining accurately its dependence on the stellar compactness. 

The paper is organised as follows. In Sect. |2]we briefly describe the initial data chosen 
and the numerical techniques employed for their evolution, while in Sect. [3] we review the 
tools used in the analysis of the data. Sect. Hcollects our results and there we first discuss the 
threshold of the bar-mode instability, its persistence, and eventually the onset of higher-mode 
dynamical instabilities. Finally, Sect. |5] contains our conclusions and the goals of our future 
research. Hereafter we use a space like signature (— , +, +, +), with Greek indices running 
from to 3, Latin indices from 1 to 3 and the standard convention for the summation over 
repeated indices. Unless explicitly stated, all the quantities are expressed in units in which 
c = G = Mq = 1. 

2. Initial data and numerical evolution method 

Our simulations involve the numerical solution in three spatial dimensions (3D) of the full set 
of Einstein equations coupled to that of a perfect-fluid matter 



and u'^ is the fluid 4-velocity, p is the fluid pressure, e the specific internal energy and 
p the rest-mass density, so that e = p{l + e) is the energy density in the rest frame of 
the fluid. The evolution of the spacetime must be supplemented by the evolution of the 
relativistic hydrodynamics equations: the conservation laws for the energy-momentum tensor 
y^jnpiy _ Q ^jjj j.jjg baryon number V^(/9ii^) = 0, complemented with an equation of state 
(EOS) of type p = p{p, e). 

The initial data for our simulations are computed as stationary equilibrium solutions for 
axisymmetric and rapidly rotating relativistic stars in polar coordinates Q. In generating 



(1) 



where 




(2) 
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Figure 1. Position on the (M/Re,f3) plane of the considered stellar models. Indicated 
respectively with stars and filled circles are the m = 2-stable and m = 2-unstable models 
belonging to the four sequences of constant rest mass. Triangles refer instead to models where 
the m=3 deformation is the fastest growing one. Indicated with a solid line is the threshold 
of the bar-mode instability, while the dashed region represents the region of our estimated 
error-bars. 



these equilibrium models the metric describing an axisymmetric relativistic star is assumed to 
have the form 

ds^ = -e^+'^dt^ + c^-'^r^ sin^ 6l(d0 - udtf + c^^{dr^ + r^dO^) , (3) 

where /x, v, uj and ^ are space-dependent metric functions. As in paper I, we assume the 
matter to be characterized by a non-uniform angular-velocity distribution of the form 

"2 r (f] - w)r2 sin^ 6'e-2'' 



i2 



1 - (f7 - w)2r2 sin^ Qq-'^^ 



(4) 



where is the coordinate equatorial stellar radius and the coefficient A is a measure of the 
degree of differential rotation, which we set to A = 1. All the equilibrium models considered 
here have been calculated using the relativistic polytropic EOS (p = with K — 100 and 
r = 2 and are members of four sequences having a constant rest mass Af* equal to 1.0 Mq, 
1.51 Mq, 2 Mq and 2.5 Mq, respectively. The main properties of the four sequences are 
reported in tables [T]|4l which report the haryonic mass Af*, the gravitational mass M, the 
angular momentum J, the rotational kinetic energy T, the gravitational binding energy W 
and the instability parameter /3 = T/|Ty|, and whose definitions are 



M : 



(-2T° - 



(5) 
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Table 1. Main properties of the stellai' models of the sequence with M, = IMq. Starting 
from the left: the name of the simulation, the compactness M/ Re, the instability parameter /3, 
the central rest-mass density pc, the ratio between the polar and the equatorial coordinate radii 
I'p/i'e, the proper equatorial radius Re, the gravitational mass M, the total angular momentum 
J divided by the square of the gravitational mass, the rotational periods at the axis Pa and at 
the equator Pe . The initial letter in the model's name indicates whether it is an unstable (U) or 
a stable (S) configuration. 



Mod 


M/Re 


/3 


pc (10-*) 


Tp/re 


Re 


M 




Pa (ms) 


Pe (ms) 


Ua8 


0.0387 


0.2814 


0.4369 


0.21044 


25.37 


0.982 


2.\42 


2.283 


4.945 


Ua7 


0.0391 


0.2790 


0.5662 


0.23177 


25.04 


0.980 


2.098 


2.207 


4.791 


Ua6 


0.0394 


0.2771 


0.6429 


0.24288 


24.83 


0.979 


2.070 


2.164 


4.702 


Ua5 


0.0398 


0.2750 


0.7087 


0.25173 


24.63 


0.979 


2.044 


2.127 


4.626 


Ua4 


0.0401 


0.2730 


0.7672 


0.25918 


24.45 


0.980 


2.021 


2.094 


4.559 


Ua3 


0.0403 


0.2711 


0.8168 


0.26550 


24.29 


0.978 


2.001 


2.068 


4.506 


Ua2 


0.0406 


0.2690 


0.8668 


0.27149 


24.11 


0.979 


1.979 


2.040 


4.449 


Ual 


0.0408 


0.2671 


0.9107 


0.27684 


23.95 


0.978 


1.960 


2.017 


4.402 


Sal 


0.0421 


0.2580 


1.0907 


0.29807 


23.22 


0.977 


1.878 


1.922 


4.209 


Sa2 


0.0424 


0.2560 


1.1270 


0.30236 


23.05 


0.977 


1.861 


1.904 


4.172 


Sa3 


0.0427 


0.2540 


1.1619 


0.30655 


22.89 


0.977 


1.844 


1.886 


4.135 



Eint = jd\^W,pe, J = jd\Tla^, (6) 

T='^jd\nTla^, W = T + Ei,,t + M,~M, (7) 

where a is the lapse function, is the square root of the three-dimensional metric 
determinant and Vt^ ~ au" is the fluid Lorentz factor. We stress that that the definitions (|5]l- 
(|7]) of quantities such as J, T, W and /3 are meaningful only in the case of stationary 
axisymmetric configurations and should therefore be treated with care once the rotational 
symmetry is lost. 

Traditionally, numerical simulations of the dynamical bar-mode instability have been 
sometimes sped up by introducing very large m=2 deformations in the initial condition. As 
discussed in paper I, the introduction of any perturbation (especially when this is not a small 
one) may lead to spurious effects and erroneous interpretations. Being aware of this, we used, 
only in some selected simulations below the threshold, initial density perturbations of the type 

Sp2{x,y,z) ^ $2 { - — P , (8) 



where 62 is the amplitude of the m=2 perturbation (which we set to be 62 — 0.01 — 0.04). 
This perturbation has then the effect of superimposing on the axially symmetric initial model 
a bar deformation that is larger than the (unavoidable) m=4-mode perturbation introduced by 
the Cartesian grid discretization. The introduction of such perturbation allowed us to estimate 
the frequency of the m=2 mode below the threshold for the onset of the instability and reduce 
considerably the computing costs in a region of the parameter space where the instability does 
not develop. 

We solve the Einstein equations ([TJ formulated as a first-order (in time) quasi-linear lH 
system of equations, where the independent variables are the three-metric 7^ and the extrinsic 
curvature A'y . In particular, we use the conformal traceless reformulation of the ADM 
system of evolution equations, first suggested in ref. ||9|, in which the evolved variables 
are the conformal factor (f>, the trace of the extrinsic curvature K the conformal 3 -metric 
jij, the conformal traceless extrinsic curvature Aij and the conformal connection functions 
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Table 2. Same quantities as in table^for the sequence of models with M, = 1.51Mq. 



Mod 


M/Re 


/3 


Pc. {U)-^) 




Rc 


M 


J/A/2 


Pa (ms) 


Pe (ms) 


Ubl3 


0.0601 


0.2812 


0.5990 


0.20012 


24.31 


1.462 


1.753 


1.723 


3.910 


Ubl2 


0.0622 


0.2761 


0.9938 


0.24151 


23.52 


1.462 


1.679 


1.599 


3.655 


Ubll 


0.0626 


0.2743 


1.0920 


0.25012 


23.31 


1.460 


1.660 


1.572 


3.598 


UblO 


0.0633 


0.2721 


1.1960 


0.25858 


23.08 


1.461 


1.639 


1.542 


3.536 


Ub9 


0.0638 


0.2701 


1.2844 


0.26554 


22.88 


1.460 


1.621 


1.517 


3.486 


Ub8 


0.0642 


0.2686 


1.3465 


0.27028 


22.73 


1.460 


1.608 


1.500 


3.451 


Ub7 


0.0646 


0.2671 


1.4055 


0.27474 


22.59 


1.459 


1.596 


1.485 


3.418 


Ub6 


0.0651 


0.2651 


1.4812 


0.28033 


22.40 


1.459 


1.579 


1.465 


3.377 


Ub5 


0.0656 


0.2631 


1.5534 


0.28560 


22.22 


1.459 


1.564 


1.446 


3.339 


Ub4 


0.0659 


0.2621 


1.5879 


0.28813 


22.13 


1.458 


1.557 


1.437 


3.321 


Ub3 


0.0664 


0.2595 


1.6730 


0.29433 


21.91 


1.456 


1.539 


1.416 


3.278 


Ub2 


0.0669 


0.2581 


1.7233 


0.29779 


21.78 


1.457 


1.527 


1.403 


3.251 


Ubl 


0.0674 


0.2551 


1.8120 


0.30450 


21.54 


1.452 


1.509 


1.384 


3.210 


Sbl 


0.0682 


0.2541 


1.8600 


0.30691 


21.42 


1.461 


1.497 


1.368 


3.179 


Sb2 


0.0682 


0.2530 


1.8845 


0.30915 


21.35 


1.456 


1.492 


1.364 


3.171 


Sb3 


0.0684 


0.2520 


1.9155 


0.31134 


21.27 


1.456 


1.485 


1.357 


3.156 


Sb4 


0.0687 


0.2503 


1.9620 


0.31500 


21.14 


1.452 


1.476 


1.348 


3.137 


Sb5 


0.0703 


0.2451 


2.1280 


0.32600 


20.70 


1.456 


1.439 


1.308 


3.057 


Sb6 


0.0713 


0.2403 


2.2610 


0.33600 


20.32 


1.449 


1.411 


1.282 


3.002 



Table 3. Same quantities as in table[T]for the sequence of models with A/, = 2Mq. 



Mod 


M/Re 


13 


Pc (10-*) 




Rc 


M 




Pa (ms) 


Pe (ms) 


Uc6 


0.0841 


0.2790 


0.9669 


0.21142 


22.82 


1.920 


1.495 


1.317 


3.161 


Uc5 


0.0858 


0.2761 


1.2663 


0.23377 


22.32 


1.916 


1.460 


1.260 


3.042 


Uc4 


0.0880 


0.2716 


1.6079 


0.25516 


21.74 


1.913 


1.420 


1.200 


2.916 


Uc3 


0.0900 


0.2670 


1.8982 


0.27117 


21.23 


1.911 


1.384 


1.152 


2.815 


Uc2 


0.0917 


0.2631 


2.1264 


0.28298 


20.82 


1.908 


1.356 


1.116 


2.741 


Ucl 


0.0931 


0.2595 


2.3176 


0.29240 


20.48 


1.907 


1.333 


1.088 


2.682 


Scl 


0.0970 


0.2500 


2.8043 


0.31526 


19.61 


1.902 


1.275 


1.021 


2.543 


Sc2 


0.0980 


0.2480 


2.9091 


0.31987 


19.42 


1.902 


1.263 


1.007 


2.514 



r*. The solution of the hydrodynamics equations is obtained by using the general-relativistic 
hydrodynamics code Whisky, in which the hydrodynamics equations are written as finite 
differences on a Cartesian grid and solved using high-resolution shock-capturing schemes, as 
described in ref. ITOl . During the evolution we use the "ideal-fluid" EOS: p ~ {T—l)pe. Full 
details of the numerical scheme and the gauge conditions used are reported in paper I. 



3. Methodology of the analysis 

A number of different quantities are calculated during the evolution to monitor the dynamics 
of the instability. Among them is the quadrupole moment of the matter distribution, which we 
compute in terms of the conserved density -^M^ p rather than of the rest-mass density p or 
of the Too component of the stress energy momentum tensor 

P'' = j (fx ^W^px^x^ . (9) 

Of course, the use of -^7!^ p in place of p or of Too is arbitrary and all the three expressions 
would have the same Newtonian limit. However, we prefer the form (|9]l because y^Wj^ p 
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Table 4. Same quantities as in table^for the sequence of models with M, = 2.5Mq. 



Mod 


M/Re 


/3 


Pc (lO-*) 




Re 


M 


J/A/2 


Pa (ms) 


Pe (ms) 


Udl5 


0.1113 


0.2771 


1.3117 


0.21259 


21.29 


2.369 


1.320 


1.025 


2.630 


Udl4 


0.1124 


0.2760 


1 .4608 


0.22121 


21.07 


2.368 


1.307 


1.005 


2.588 


Udl3 


0.1132 


0.2750 


1.5826 


0.22785 


20.90 


2.366 


1.298 


0.990 


2.555 


Udl2 


0.1140 


0.2740 


1.6961 


0.23372 


20.74 


2.364 


1.289 


0.976 


2.526 


Udll 


0.1149 


0.2730 


1.8103 


0.23925 


20.58 


2.364 


1.280 


0.962 


2.496 


UdlO 


0.1156 


0.2721 


1.9104 


0.24399 


20.44 


2.362 


1.272 


0.950 


2.471 


Ud9 


0.1164 


0.2710 


2.0153 


0.24868 


20.29 


2.362 


1.264 


0.938 


2.445 


Ud8 


0.1171 


0.2701 


2.1102 


0.25286 


20.16 


2.360 


1.256 


0.927 


2.423 


Ud7 


0.1186 


0.2680 


2.3022 


0.26083 


19.89 


2.359 


1.241 


0.906 


2.377 


Ud6 


0.1189 


0.2675 


2.3534 


0.26291 


19.82 


2.358 


1.238 


0.901 


2.366 


Ud5 


0.1196 


0.2666 


2.4336 


0.26609 


19.71 


2.357 


1.231 


0.892 


2.348 


Ud4 


0.1207 


0.2650 


2.5698 


0.27128 


19.52 


2.356 


1.221 


0.878 


2.318 


Ud3 


0.1220 


0.2631 


2.7402 


0.27760 


19.29 


2.353 


1.208 


0.861 


2.282 


Ud2 


0.1254 


0.2580 


3.1583 


0.29210 


18.74 


2.349 


1.178 


0.821 


2.198 


Udl 


0.1302 


0.2510 


3.7335 


0.31027 


18.00 


2.343 


1.138 


0.772 


2.094 


Sdl 


0.1314 


0.2491 


3.8899 


0.31503 


17.81 


2.341 


1.127 


0.760 


2.068 


Sd2 


0.1321 


0.2480 


3.9735 


0.31745 


17.71 


2.340 


1.122 


0.753 


2.054 


Sd3 


0.1329 


0.2470 


4.0607 


0.31991 


17.61 


2.340 


1.116 


0.746 


2.040 


Sd4 


0.1341 


0.2450 


4.2188 


0.32450 


17.42 


2.337 


1.106 


0.735 


2.016 


Sd5 


0.1376 


0.2400 


4.6380 


0.33596 


16.95 


2.333 


1.081 


0.706 


1.955 


Sd6 


0.1412 


0.2350 


5.0710 


0.34714 


16.49 


2.329 


1.056 


0.678 


1.897 



is a quantity whose conservation is guaranteed by the form chosen for the hydrodynamics 
equations. The quantity (|9|l can be conveniently used to quantify both the growth time of the 
instability T2 and the oscillation frequency of the unstable bar-mode once the instability is 
fully developed /2. (Hereafter we will indicate respectively with Ti and fi the growth time 
and frequency of the m = i unstable mode.) 

In practice, we perform a nonlinear least-square fit of the computed quadrupole I^^{t) 
and we generally use the xy component, with the trial functions 

P\t) = (/^■^)o e*/^^ cos(27r/2t + M . (10) 

Furthermore, we define the modulus I{t) of the two components of the quadrupole in the xy 
plane and the distortion parameter 77 (t) as 

/ = ]-J(2i^vY + (i^^ - ivvy- , 77=——^——. (11) 

and the instantaneous orientation of the bar is given by 

ixx _ lyy 



bar = tan-1 ( I . (12) 



Finally, as a useful tool to describe the nonlinear properties of the development and 
saturation of the instability, the rest-mass density is decomposed into its Fourier modes Pm {t)'- 



P,n= Id^xpe'""^ . (13) 



The phase = arg(P,„) essentially provides the instantaneous orientation of the 7T! -th mode 
when the corresponding mode has a nonzero power. Note that despite their denomination, the 
Fourier modes (fTSI l do not represent proper eigenmodes of oscillation of the star While, in 
fact, the latter are well defined only within a perturbative regime, the former simply represent a 
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Table S. Least-squai'e fit of the value of f) at the threshold for the development of the bar-mode 
instability for the four series of models reported in table [8] The critical value for the onset of 
the instability /3c is the value of (3 for 1/t| = (r is measured in ms) and the digits in 
brackets represent the en'or in the fit. The results of the fits ai'e shown in figure[3] 



M, 


= 1.0 Mq 


p 


= 0.2598(8) 


+0.0379(19) 


(1/^2)2 


M* 


= 1.5 Mq 


p 


= 0.2558(5) 


+0.0236(8) 




M* 


= 2.0 Mq 


p 


= 0.2528(15) 


+0.0161(13) 


(1/^2)2 


M. 


= 2.5 Mq 


p 


= 0.2494(14) 


+0.0116(8) 


(1/^2)2 



Table 6. Least-square fit of the value of the frequency /b (in Hz) of the bar-mode (Equation 
not ) at the threshold for the onset of the bai-mode instability as a function of 9 = (/S — 13c)/ Pc 
for the four series of models at constant baryonic mass. The value of the frequency /2 at the 
threshold is the value for 9 = and the digits in brackets represent the error in the fit. 





= 1.0 M0 


/b 


= 384(7) 


-8(25) 10 6» 


-6(2) 1000 02 


M. 


= 1.5 Mq 


/b 


= 551(8) 


-4(3) 100 


-6(2) 1000 02 


M* 


= 2.0 Mq 


/b 


= 738(32) 


-7(11) 100 


-6(8) 1000 02 


M* 


= 2.5 Mq 


/b 


= 991(10) 


-10(3) 100 


-9(2) 1000 02 



tool to quantify, within the fully nonlinear regime, what are the main components of the rest- 
mass distribution. Stated differently, we do not expect quasi-normal modes of oscillations 
to be present but in the initial and final stages of the instability, for which a perturbative 
description is adequate. 

While all quantities (l9])-(fT3]) are expressed in terms of the coordinate time t and do not 
represent therefore invariant measurements, the lengthscale of variation of the lapse function 
at any given time is always larger than twice the stellar radius at that time, ensuring that events 
on the same timeslice are also close in proper time. As representative examples, we note that 
for the most compact model the values of the lapse at the centre of the star, at its surface and at 
the outer boundary are 0.67, 0.84 and 0.95, respectively. Similarly, the corresponding values 
for the least compact model are 0.92, 0.95 and 0.98, respectively. 

The simulations have been carried out on a grid with a uniform resolution of Ax/Mq ~ 
0.625 and outer boundaries at 48.75 Af0, where "radiative" boundary conditions (i.e., 
Sommerfeld outgoing boundary conditions) are applied to the field variables, while the fluid 
variables are simply not evolved. Such outer boundaries are sufficiently far from the surface 
of the star to make the use of mesh refinements in Whisky not necessary (see also Section 
VIE of paper I for a more detailed discussion of the role of the grid size on the evolution of 
the instability). Furthermore, for those models used in the extrapolation technique and that 
are largely over-critical (see Sect. 14. lb . we have imposed a "bitant" symmetry {i.e., z —z) 
and a "7r-symmetry" {i.e., a 180° rotation around the z-axis) in order to reduce the size of the 
computational domain of a factor of 4. 

4. Effects of the compactness 

4.1. Threshold of the m~2 instability 

The determination of the dependence on (3 of the frequencies and of the growth times 
of the TO = 2 bar-mode instability in the region near the threshold is particularly delicate 
as the models are only slightly over-critical, with very small growth rates and hence the 
simulations are computationally very expensive. For this purpose we here use an extrapolation 
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0.24 



0.26 



0.28 



Figure 2. Critical diagram as constructed witli tlie frequencies and growth times relative to the 
unperturbed models of table[8] The sohd lines represent the two fitted curves for r2(/3) and 
T2{f3), while the dotted Hues the corresponding extrapolations below the threshold. Triangles 
refer to the unperturbed models of table [8] that have not been used for the fit, squares to the 
perturbed models of table |7] (squares), and open circles to the the models dominated by the 
m=3 deformation. 



technique already described in paper I, where it was shown to be both accurate and 
robust. In essence, we exploit the results of the classical Newtonian study of the bar-mode 
instability of Maclaurin spheroids of incompressible and self-gravitating Newtonian fluid 
in equilibrium ifTTl . extrapolating, via suitable fits, its predictions to a general-relativistic 
context. We recall, in fact, that in the classical scenario the eigenfrequency of the m=2 bar- 
mode can be expressed in terms of two real and differentiable functions of /3; ft = 27r/2 and 
1 /r|, in a relation of the type 

The bar-mode becomes unstable when the function 1/t| changes sign, with the square root 
going from being imaginary to being real. The value of f3 at which this change of sign happens 
represents then the threshold for the onset of the instability (3c; clearly, for models above the 
threshold, H/Att and T2 represent the pattern speed and the growth time of the unstable bar 
deformation of the considered star model, respectively. 

With the rather reasonable assumption that the two functions 51 and 1 /r| are regular 
also in full General Relativity, we expand them in a Taylor series around the threshold and 
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Table 7. For models below the threshold for the onset of the m=2 bar-mode instability the 
table reports the measured frequencies J2 of the bar deformation as well as the frequencies fz 
and growth times rs of the m=3 deformation. For each model we also indicate the value of /3, 
the grid resolution Ax/Mq and the initial bar-mode perturbation 52- All models have been 
evolved without a 7r-symmetry. 



Model 


/3 


Ax/Mq 


S2 


/2 (Hz) 


/3 (Hz) 


T-A (ms) 


Sa3 


0.254 


0.5 


0.30 


424+5!^ 
—62 






Sa2 


0.256 


0.5 


0.01 


409+ 5n 

40liS 






Sal 


0.258 


0.5 


0.01 




* 


Sb6 


0.240 


0.5 


0.04 


736+!S 






Sb5 


0.245 


0.5 


0.04 


694+1^ 

50 






Sb4 


0.250 


0.5 


0.04 


636+^!^ 






Sb3 


0.252 


0.5 


0.04 








Sb2 


0.253 


0.5 


0.04 


58811 






Sbl 


0.254 


0.5 


0.04 


57611 






Ubl 


0.255 


0.5 


0.04 


56411 






Sbl 


0.254 


0.5 


0.0 


6O6I21 


1302;ii 


3.8 


Ubl 


0.255 


0.5 


0.0 


586li2 


1284+12 


10.0 


Scl 


0.250 


0.625 


0.01 


782+^0 

' " — 70 






Sc2 


0.248 


0.625 


0.01 


818+1° 






Sc2 


0.248 


0.625 


0.0 




1746^^? 


4.0 


Sd6 
Sd4 


0.235 

U.Z4U 

0.245 


0.625 

U.Dzj 

0.625 


0.01 

U.Ui 

0.01 


1400I24 

i3U4_gg 

126211 


245411 

OQI 0+122 

224811 


2.7 

1 A 

2.9 


Sd3 
Sdl 


0.247 
0.249 


0.625 
0.625 


0.01 
0.01 


109811 
104811 


22OOII 
21361^4 


2.5 

3.5 


Sd6 


0.235 


0.625 


0.0 


* 


2474l^« 


2.7 


Sd5 


0.240 


0.625 


0.0 


* 


2372II 


2.7 


Sd4 


0.245 


0.625 


0.0 


* 


22961^^ 


2.2 


Sd2 


0.248 


0.625 


0.0 




222811 


2.5 


Sdl 
Udl 


0.249 
0.251 


0.625 
0.625 


0.0 
0.0 


144611 
IO22II 


217411 
215811 


2.2 
2.9 



* A dynamical m=3 deformation instability cannot be detected. 

* The frequency of the m=2 bar-mode cannot be measured. 



expressed them in terms of five unknown coefficients /c, /i^"* , /i^"* , fc, Pc, i e. 



(15) 



\~^{P-(ic) + 0{{l3-(i,f). (16) 

Expressions ( fTSI ) and ( fTSt represent very good approximations to the actual data and 
the five parameters can be determined straightforwardly by fitting the pattern speeds and the 
growth times obtained in the largely over-critical models. For these models, we recall, the 
development of the 111=2 bar-mode deformation is very rapid, the extraction of the instability 
parameter is robust and it can be safely simulated even at rather low resolutions (see paper 1). 

In practice, using the data obtained from the simulations of the four sequences of initial 
models with constant baryonic mass, we have computed the values for /2 and T2 by means of 
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Table 8. For models above the threshold the table reports the measured frequencies /2 and 
growth times T2 of the m = 2 deformation as well as the value of /3. All models have been 
evolved with a grid resolution Ax/AIq = 0.625 and no initial perturbation; furthermore a 
TT-symmetiy has been used for all models with the exception of models Ubl and Udl that are 
closer to the threshold. 



Model 


/3 


T2 (ms) 


/2 (Hz) 


Model 


/3 


T2 (ms) 


/2 (Hz) 


Ual 


0.2671 




377I4 


Ucl 


0.2596 


1.564_o 17 


■71 v+J-0 
717„g 


Ua2 


0.2690 


2-019_o.i3 


3731'^ 


Uc2 


0.2631 


-, 0-77+0.07 


694lg* 


Ua3 


0.27 11 


1 CQQ + 0.11 


36813 


Uc3 


0.2670 


I.O55I0 10 


679I7 


Ua4 
Ua5 


0.2730 
0.2750 


1.702lg:J0 


36411 
359I4 


Uc4 
Uc5 


0.2716 
0.2761 


0.915l°:°«2 
0.8291°:°^ 


650I5 
6211^ 


Ua6 


0.2771 


l-482_o 


352I4 


Uc6 


U.I /yu 


0-'95_o 03 




Ua7 


0.2790 


1 /inc+O.ll 
l-405_o(,4 


345 _ 2 


Udl 


0.2510 


^ 4.2 


1009_ye 


Ua8 


0.2814 


-, qi O+O.IO 


335I2 


Ud2 


0.2580 


1.193I0 11 


9451 j^4 


Ubl 


0.2551 




58d1j^2 


Ud3 


0.2631 


0.938Io(,2 


9071^ 


UDZ 


U.zjoi 


r,, C- + 2.62 
O.Z10_Q 


c:c:9+40 


Ua4 


U.zOJU 


0.846_o 03 


SQ1 +4 

8yi_4 


Ub3 


0.2595 


2.758 


5441^:^ 


Ud5 


0.2666 


n onn+0.07 

0-800-0.05 


88lly 


Ub4 


0.2621 


1.920i°:^^ 


540l| 


Ud6 


0.2675 


0.8131°°^ 


87ll^^ 


Ub5 
Ub6 
Ub7 


0.2631 
0.2651 
0.2671 


1.803in^ 
1 441 +010 


5351^ 
5271^ 
520I* 


Ud7 
Ud8 
Ud9 


0.2680 
0.2701 
0.2710 


n 774+0.07 
n 7oi;+0.05 
794+0.05 


865I4 

8471L 

8321^ 


Ub8 


0.2686 


1.000_g Q3 


51611 


UdlO 


0.2721 


0.7161°:°^ 


826I7 


Ub9 


0.2701 


1 2qo+''°'' 

l.ZaU_Q Qg 


509l| 


Udll 


0.2730 


0.6971°:°^ 


8I6I9 


UblO 


0.2721 


1 214+°°* 


502I4 


Udl2 


0.2741 


0.690l°:°3 


804+ g 


Ubll 


0.2743 


1.1261°:°° 


492I5 


Udl 3 


0.2750 


0.679l°:°5 


793l*i 


Ubl2 


0.2761 


1.0871°:°! 


483+* 


Udl4 


0.2760 


0.660+°°® 


782lg 


Ubl3 


0.2812 


0.959l°:J? 


4501^ 


Udl 5 


0.2771 


0.6541°;°^ 


767lg 



* The growth rate of the m=2 bar-mode cannot be reliably measured for this model. 



a nonlinear least-square fit to the trial form of eq. ( fTOb . Making use of these results, which are 
collected in table[8] we have then computed the unknown coefficients fc, fc^\ f'^\ k and /3c 
through a least-square fit of the values for /2 and T2 once expressed as functions of (3 as in 
eqs. ( flSl ) and ( fT6] l. 

The results of these fits are reported in Tables [5] and |6] and summarised in Fig. |2] In 
particular, for each of the four sequences this figure shows with solid lines the two fitted curves 
for n(/3) and T2(/3), and with dotted lines the corresponding extrapolations for models below 
the threshold. In addition, different symbols are used to mark the results of the numerical 
simulations, with "bare" error bars denoting the unperturbed unstable models (as reported in 
Table [8] and which have been used for the fits), triangles denoting the unperturbed models (as 
reported in Table [8] and not used for the fits because of the large error in determining their 
evolution parameters), squares denoting the stable perturbed models (as reported in Table |7]) 
and open circles denoting the models dominated by the m=3 instability (and again reported 
in Table |7]). It is worth noting that for these last models the frequency of the m=2 mode is 
considerably altered by the growing m = 3 deformation, which rapidly dominates the small 
initial m=2 bar-mode perturbation used for the simulation of these models. 

We also note that the error bars in Table [8] are computed in different ways for the 
growth rates and for the frequencies. In the first case they are computed as the difference 
between the minimum and maximum values of dlog{I{t))/dt in the time intervals in which 
the quadrupolar deformation lies between the 5% and the 35% of its maximum amplitude. 
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Figure 3. Left panel: Indicated with triangles are the data already reported in Fig. [2] but 
shown here as a function of l/r^ for the four constant-rest-mass sequences considered. The 
open circles indicate instead the extrapolation to Pc- The numerical values of the four fits are 
reported in table[5] Right panel: The open circles are the same as in the left panel but shown as 
a function of the baryonic mass Mt ; the filled circle, on the other hand, represents the hmit of 
zero rest-mass and can be compared with the threshold value calculated in Newtonian gravity 
and with a linear stability analysis in ref. 1131 . 



In the second case, instead, the error bars are determined using the minimum and maximum 
values in the time interval over which the pattern speeds are extracted from the collective 
phase (/)bar(i) [c/ eq. (fT2]il. The frequencies and the error bars for the stable models reported 
in Table [T] and Fig. |2l on the other hand, are computed using the Lomb's power spectrum 
analysis lfT2ll . which is better suited to study a signal spanning over a short time interval and 
which is comparable with the main frequency of the Fourier transform. In particular, they 
refer to the frequency at which the Fourier transform of the imaginary part of P,„ (t) has its 
maximum, while the error bars to the interval in frequencies where the Fourier transform is 
above 1/2 of its maximum value. 

The left panel of Fig. [3] offers a different view of the data already reported in Fig. |2] 
by showing /3 as a function of 1 /r^ for the four sequences considered and thus highlighting 
the very good approximation in the ansatz ( fT6l l. Indicated with open circles are also the 
extrapolations of the critical value /3c and these are then reported in the right panel of Fig. [3] 
as a function of the baryonic mass Af* of the four sequences. Clearly, they show a linear 
dependence on AI^,, which can be expressed phenomenologically as 



and allows us to extrapolate the value of the threshold for in the limit of a zero baryonic mass. 
Interestingly, the resulting number (3c{M^, = 0) = 0.266 is in very good agreement with the 
value of /3c = 0.266 obtained in Newtonian gravity ifTSll but through a linear stability analysis 
for a sequence of equilibrium models with the same polytropic index and degree of differential 
rotation used here. This agreement represents an additional confirmation of the accuracy and 
robustness of our extrapolation method in determining the position of the threshold. 

Using the phenomenological dependence of the threshold (3c on the stellar rest mass 
given by expression ( fTTI i. we have also reconstructed the dependence of /3c on the stellar 
compactness M/Re- In practice, for a large number of values of A/* between and 2.5Mq we 




(17) 
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Figure 4. Comparison of the evolution of the distortion parameters r} for a low-mass, low- 
compactness model Ub4 and a high-mass, high-compactness model Ud2. Indicated with a 
dotted and soHd line are the simulations for Ub4 and Ud2 using a 7r-symmetry, while indicated 
with a dashed lines is the evolution for model Ud2 without 7r-symmetry. All the evolutions 
have been suitably shifted in time so as to have the maximum deformation at t = 0. 

have computed the value of the proper equatorial radius i?e and of the gravitational mass M 
of the corresponding stellar model in equilibrium and then used eq. ( fTTI ) to estimate the value 
of /3c- The result of this is shown as a solid line in Fig.[Tl with the dashed band representing 
the estimated error obtained using the least-square fitting. Furthermore, a good polynomial 
reconstruction of the median of the error bar in Fig. [T] suggests a quadratic dependence of the 
threshold on the compactness, with coefficients given by 

/3c = 0.266 -0.18 (^1^^ +0.36 (^^^ . (18) 

Once perturbative calculations will be developed in the regime of rapid and differential 
rotation considered here, expression (fTsT i can be used as a guideline to the perturbative 
approach and the numerical measurements of the threshold used to assess the validity and 
accuracy of the perturbative approximation. 

Overall, these results represent the first quantitative determination of the dependence of 
the critical jS on the compactness of the star for a selected profile of differential rotation and 
EOS. A similar investigation was carried out also in ref. lfT4ll using the post-Newtonian (PN) 
approximation and rather large initial m=2 perturbations with 52=0.1. While the results of 
that analysis had the limitation of the PN approximation and were affected by the large initial 
perturbations, they also provided the first evidence that the threshold is smaller for stellar 
models with larger compactness. 
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4.2. Persistence of the bar 

As mentioned in the Introduction and for all of the models considered here, any non- 
axisymmetric deformation that develops as a result of a dynamical instability is also 
suppressed over a time comparable {i.e., of the order of a few) dynamical timescales. In 
addition, we have found that the persistence of the bar deformation increases as /3 — /3c tends 
to zero. Besides confirming what already discussed in paper I, this behaviour matches the 
expectation that the persistence of the bar is related to the degree of overcriticality, with the 
duration of the saturation increasing as the threshold is approached. For a star with a small 
compactness this time would tend to the radiation-reaction timescale for a model with (3 = pc 
and to zero for a model with /3 ^ Pc- For a star with large compactness, on the other hand, 
the persistence near the threshold can be further reduced by the stronger gravitational fields. 

Here, we also intend to gain insight on the role played by the stellar compactness on the 
persistence of the bar deformations and, to this scope, we have extended the time over which 
the simulations are carried out for some selected models. Figure |4] summarizes the results 
of these extended simulations by reporting the evolution of the distortion parameter -q for two 
models, i.e., Ub4 (dotted line) and Ud2 (solid line) having different masses and compactnesses 
{i.e., M = 1.5 Mq, M = 2.5 Mq and M/R,, ~ 0.066, M/Re 0.125, respectively), but 
with a similar distance from the threshold {i.e., (3 — Pc — 0.007 — 0.008). The two evolutions 
have been suitably shifted in time so as to have the maximum deformation at t = 0. 

Both of these models have been evolved using a 7r-symmetry in order to remove the 
effect of the odd-mode coupling on the suppression of the bar deformation and yet they show 
a remarkable difference. The low-mass, low-compactness model Ub4 reaches and maintains 
a rather large bar deformation over several milliseconds of evolution. The high-mass, high- 
compactness model Ud2, on the other hand, reaches comparable deformations but these are 
rapidly suppressed over a few milliseconds, despite the use of the 7r-symmetry. This behaviour 
underlines something that was already remarked in paper I, namely that, depending on the 
specific stellar properties, the bar deformation can be suppressed also by factors other than 
the mode coupling which is most effective for models near the threshold. In particular, for 
low-mass models (such as model U13 in paper I) this effect is simply the excess of rotational 
kinetic which is efficiently converted into internal one (see discussion in Sect. VIB of paper I). 
For high-mass and high-compactness models (such as model Ud2 above), on the other hand, 
the strength of the gravitational field, together with the excess kinetic energy, are very efficient 
in suppressing the bar-deformation. This was not evident in the stellar models considered in 
paper I, which all had A/* = 1.5 Mq. 

To confirm that the m = 1-mode coupling plays no significant role in the dynamics of 
model Ud2, we show with a dashed line in Figure |4]the evolution of the bar deformation for 
model Ud2 in a simulation in which the 7r-symmetry was not enforced. As it is obvious in 
the comparison between the solid and dashed lines, the lack of symmetry does not change the 
suppression of the bar-mode deformation and a difference emerges only at i 5 ms, when the 
simulation with 7r-symmetry tends to revive the bar deformation, while this does not happen 
in the simulation without the 7r-symmetry. 

4.3. Unstable deformations with m=3 

While studying the dynamics of models with > 2Mq and values of P near the threshold 
for the development of the bar-mode instability, we have also found stellar models that show 
the development of a dynamical instability with a dominant m=3 deformation. Interestingly, 
this instability developed without the introducing of any m = 3 initial perturbation as was 



Dynamical non-axisymmetric instabilities in rotating relativistic stars 



14 



; t = 20.8 ms 

r 


- 


^ f 






80% /)/ /,-■■',■■:';■,-■■ i':: I 







10 20 30 




Figure 5. Left panel: Rest-mass isodensity model Sd4. Starting from the center the sohd hnes 
represent respectively: 80%, 70%, 60% and 50% of the maximum density. Dotted lines are 
1/2" times the maximum density, with n = 2, . . . , 11. The snapshot was taken at t=20.8 ms. 
Rigth panel: Evolution of the global rest-mass density modes defined in eq. 1131 for the same 
simulation. 



instead done in ref. ifTSl . More precisely, we found four models (Sc2, Sd6, Sd5, Sd4) for 
which the m~?> mode is the fastest growing dynamically unstable deformation. In addition, 
we also see evidence of a dynamical unstable growth of the to=3 mode for models Sbl, Ubl, 
Sd3, Sd2, Sdl and Udl. However, because all the above models lay very close to the threshold 
for the onset of the bar-mode instability, nonlinear mode-couplings may be very important in 
this region and the cause of the observed growths. 

The general behavior of the mode dynamics is reported in Fig.|5]for the representative 
model Sd4, where in the left panel we show the rest-mass isodensity contours in the equatorial 
plane and at t = 20.8 ms, while in the right one the evolution of the global modes as defined 
in eq. ( fTSl ) but as measured with respect to the centre of mass of the system (see discussion in 
Sect. IV of paper 1 for the definition of the centre of mass in this context). When looking at this 
panel it is evident that the m=3 unstable deformation is the fastest growing one. Furthermore, 
the m=2 does not show any sign of an unstable growth, which is instead shown by the m=l 
deformation. Such a growth starts roughly at the time when the m=l and m=2 modes have 
powers comparable with the background m = 4, i.e. at t ~ 20 ms, and it continues for the 
rest of the simulation. The results on the frequencies /s and the growth times T3 of the m=3 
instability are collected in Table|7]and have been computed using the same analysis discussed 
for the m=2 deformation. It should be noted that the properties of the m=3 instability are 
much more difficult to estimate accurately as the corresponding deformations are smaller than 
those measured for the bar-mode instability. 

It is probably worth stressing that this is the first time that similar instabilities are seen in 
fully general relativistic simulations of stellar models with a stiff EOS and a moderate degree 
of differential rotation. In a recent work carried out in Newtonian physics [|6l, in fact, the 
development of deformations with m = 3 was also found but for stellar models with lower 
values of j3 and with a stronger differential rotation. Furthermore, because the models were 
studied by imposing an initial artificial m=2 perturbation that could alter the global dynamics, 
it is indeed difficult to understand if the origin of this unstable deformation is purely dynamical 
or triggered by a nonlinear coupling of modes. Nevertheless, it is interesting to note that in 
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all the simulations reported in JS) as well as in those presented here the m=l deformation 
becomes the dominant one in the final stages of the evolution. This lends support to the 
idea that the suppression of non-axisymmetric deformations over a dynamical timescale is a 
generic feature of these instabilities for isolated stars and not necessarily restricted to stellar 
models with high values of (3. 

At the moment and besides the classical m=2 bar-mode instability, a proper 
understanding of the conditions that lead to the development of dynamical instabilities is 
still lacking. The first attempts of interpreting these instabilities, and especially the low-/3 
m=l instability, have been made recently by several authors ||6][T6|. However, further work 
is needed to clarify the role of j3, of the differential rotation law, of the EOS and of the 
compactness in determining the growth times, the maximum amplitudes and the persistence 
of these unstable deformations. Most importantly, work will be needed to finally reach 
sufficient conditions for the onset of these instabilities whose development has been revealed 
by numerical simulations. 

5. Conclusions 

We have presented accurate simulations of the bar-mode dynamical instability in full General 
Relativity. The main motivation behind this work was to address some important open 
questions about the nonlinear features of non-axisymmetric dynamical instabilities in rapidly 
rotating compact stars. The most important among these questions, because of the impact 
it has on the global detectability of these stars as sources of gravitational waves, is the 
determination of the timescale over which the non-axisymmetric deformations persist once 
these are fully developed. 

In order to reach a better understanding of the physics governing dynamical instabilities, 
we have analysed the onset and development of the bar-mode instability for a large number 
of stellar models spanning a wide range of masses and angular momenta. The initial models 
have been calculated as stationary equilibrium solutions for axisymmetric and rapidly rotating 
neutron stars modeled as relativistic polytropes with adiabatic index F = 2 and polytropic 
constant K = 100. All the stars have been constructed with a differential-rotation profile 
having A = 1 and as members of four sequences of constant rest-mass, with A/, = IA/0, 
I.SIM©, 2Mq and 2.5Mq, respectively. This large set of initial data containing a total of 59 
models has allowed us to confirm and extend the results presented in paper 1. More precisely, 
we have analysed the effects that the stellar compactness has on the threshold for the onset of 
the dynamical bar-mode instability as well as of other dynamical instabilities. Moreover, using 
an extrapolation technique developed and tested in paper 1, we have determined the threshold 
with great accuracy and for a wide range of compactnesses, finding a good agreement with 
the Newtonian prediction and improving the previous PN estimates made in ref. lfT4l . 

While studying the dynamics of models with > 2Mq and values of (3 near the 
threshold for the development of the bar-mode instability, we have also found stellar models 
that show the development of a dynamical instability with a dominant m = 3 deformation, 
and without any seed perturbation of that type. The appearance of these instabilities, whose 
growth-time and frequency have been computed using the same methodology developed for 
the bar-mode instability, may be rather generic in stars with high mass and deserve additional 
attention. 

Finally, we remark that for all the simulated models, the deformations generated by the 
non-axisymmetric instabilities are suppressed over a dynamical time-scale either as a result 
of nonlinear mode-couplings or as a result of the conversion of the excess rotational kinetic 
energy into internal one. In all cases we have observed the emergence of a residual m=l 
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deformation in the final stages of the instability and before an axisymmetric configuration 
is recovered. These results confirm our previous findings presented in paper 1 and are in 
agreement with those recently reported in |[6il, thus lending support to the idea that the 
suppression of non-axisymmetric deformations over a dynamical timescale is a generic feature 
in isolated stars. Further work is clearly needed to confirm or reject this conjecture and to 
derive the sufficient conditions for the onset of these "odd-m" instabilities. 

Acknowledgments 

It is a pleasure to thank Harald Dimmelmeier and Ian Hawke for useful discussions and 
suggestions. GMM is also grateful to Jan Christian Bryne and the Bergen Center for 
Computational Science for the kind hospitality. Support for this research comes also through 
the SFB-TR7 of the German DFG and through the 0G51 of the Italian INFN. GMM is 
presently a "Dalla Riccia" fellow at the MPA. All the computations were performed on the 
cluster for numerical relativity "Albert" at the University of Parma and "Peyote" at the AEI. 

References 

[1] Baiotti L, De Pietri R, Manca G M and Rezzolla L 2007 Pliys. Rev. D 75, 044023. 
[2] Shibata M and Sekiguchi Y 2005 Phys. Rev. D 71, 024014. 
Saijo M 2005, PIm. Rev D71, 104038. 

Ott C D, Ou S, Tohline J E and Burrows A 2005 Astrophys. J. 625, LI 19-L122. 

Ott C. D. et al 2006 "3D CoUap.se of Rotating Stellar Iron Core.s in General Relativity with Microphysics" , 
ai-Xiv:astro-plV0609819 A. 
[3] Dessart L, Burrows A, Ott C D, Livne E, Yoon S C and Langer N 2006 Astrophys. J. 644, 1063. 
[4] Shibata M, Taniguchi K and Uryu K 2003 Phys. Rev D 68, 084020. 

Shibata M, Taniguchi K and Uryu K 2005 Phys. Rev D 71, 084021. 

Shibata M and Taniguchi K 2006 Phys. Rev D 73, 064027. 
[5] Bonaldi M et al 2003 Phys. Rev D 68, 102004. 

Bonaldi M et al 2006 Phys. Rev D 74, 022003. 
[6] Ou S and TohHne J E 2006 Astrophys. J. 651, 1068. 
[7] Stergioulas N and Friedman J L 1995 Astrophys. J. 444, 306. 

[8] Richtmyer R D and Morton K W 1967 Difference Methods for Initial Value Problems. Interscience Publishers, 
New York, 1967. 

[9] Nakamura T, Oohara K and Kojima Y 1987 Prog. Theor Phys. SiippL, 90, 1. 
[10] Baiotti L, Hawke I, Montero P and Rezzolla L 2003 A new three-dimensional general-relativistic 
hydrodynamics code, in R. Capuzzo-Dolcetta, editor. Computational Astrophysics in Italy: Methods and 
Tools, volume 1, page 327, Trieste, 2003. Mem. Soc. Astron. It. Suppl. 
Baiotti L et al 2005 Phys. Rev D 71, 024035.. 
[11] Chandrasekhar S 1969 Ellipsoidal Figures of Equilibrium. Yale Univ. Press, New Haven. 
[12] Lomb N R 1976 A.str Space Set, 39:447^62, 1976. 
[13] Karino Sand Eiiguchi Y 2003 Astrophys. J. 592, 1 1 19-1 123., 2003. 
[14] Saijo M, Shibata M, Baumgarte T W and Shapiro S L 2001 Astrophys. J. 548, 919-931. 
[15] Rampp M, Mueller E and Ruffert M 1998 Astron. Astrophysics, 332, 969. 
[16] Watts A L, Andersson N, and Jones D I 2005 Astrophys. J. 618, L37. 

Saijo M and Yoshida S 2006 Mon. Not. Roy Astron. Soc. 368, 1429-1442. 



